library(raster)

library(plyr)
library(dplyr)
library(reshape)

library(ggplot2)
library(ggpmisc)
library(ggpubr)

library(viridisLite)
library(colorspace)
library(RColorBrewer)

library(sf)
library(stringr)

memory.limit(size = 160000)
## [1] Inf

#display.brewer.pal(n = 11, name ="RdYlGn")
colorvec <- brewer.pal(n = 11, name ="RdYlGn")

col.martonne <- c("Desert" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Temperate" = colorvec[8], "Humid"= colorvec[10],"Forest"= colorvec[11])
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "NA" = "#D5D8DC")

1 Data

Liens utiles:

# era5vA <- read.table("era5_AI.txt") %>%
#   mutate(tas = tas-273.15, pr = mpr * 60 * 60 * 24 * 365, rsds = ssrd / (60*60*24) * 1e-6 * 60*60*24) %>% # conversion from J/m2 to J/sm2 to MJ/m2/day  
#   select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn","rsds","ea","ET0","AI","cat.AI"))
# 
# 
# write.table(era5vA, "era5vA.txt")

era5vA <- read.table("ERA5/era5_AI.txt") %>% 
  mutate(z = z, tas = t2m - 273.15, sfcWind = wind, model = "historical", source = "ERA5") %>%  #pr, ET0 already in mm/year
  select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn","ea","z","ET0","AI","cat.AI")) %>%
  filter(lm == 1)

write.table(era5vA, "era5vA.txt")

wdcvA <- read.table("wdc_AI.txt") %>% mutate(z = z, ea = vapr, ET0 = ET0) %>% # srad in MJ/m2/day, pr and ETO already in mm/year
 select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI","cat.AI")) %>%
  filter(lm == 1)


write.table(wdcvA, "wdcvA.txt")


cmip6 <- read.table("cmip6_AI.txt")

cmip6vA <- cmip6 %>% subset(period == "1970_2000") %>%
   group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, z) %>% 
   dplyr::summarise(tas = mean(tas, na.rm = T) - 273.15, pr = mean(spr, na.rm = T), # pr in mm/year
                   sfcWind = mean(sfcWind, na.rm = T), Rn = mean(Rn, na.rm = T),
                   ea = mean(ea, na.rm = T), ET0 = mean(sET0, na.rm = T), AI = mean(AI, na.rm = T)) %>% # ET0 in mm/year
  ungroup() %>%  
  mutate(source = "CMIP6")%>%
   select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI"))


breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

cmip6vA$cat.AI <- cmip6vA$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6vA$cat.AI[which(cmip6vA$ET0 < 400)] <- "Cold"

write.table(cmip6vA, "cmip6vA.txt")
era5vA <- read.table("era5vA.txt")
wdcvA <- read.table("wdcvA.txt")
cmip6vA <- read.table("cmip6vA.txt")

2 Regions and Continents


colors_continents <- c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
                  "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf",
                  "#aec7e8", "#ffbb78", "#98df8a", "#ff9896")

ggplot(era5vA)+geom_raster(aes(x= lon, y = lat, fill = Continent))+
  scale_fill_discrete(type = colors_continents)+
  theme_void()

“Continents” to remove are: * Arctic (too little gridcells) * Atlantic (ocean) * Indian (ocean) * Pacific (ocean) * Southern (ocean)

land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
plot(land_mask)

land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land


elev <- raster("Worldclim/wc2.1_10m/wc2.1_10m_elev.tif") 
elev.df <- elev %>% projectRaster(to = land_mask) %>% as.data.frame(xy = T) %>% setNames(c("lon","lat", "z"))


ipcc_regions <- shapefile("Masks/IPCC-WGI-reference-regions-v4.shp") %>% spTransform(crs("EPSG:4326")) 

ipcc_regions.raster <- ipcc_regions %>% rasterize(land_mask)
ipcc_regions.df <- as.data.frame(ipcc_regions.raster, xy = T) %>% setNames(c("lon","lat","Continent","Type","Name","Acronym"))

ipcc_sf <- st_as_sf(ipcc_regions)

ggplot(data = filter(ipcc_sf, Type %in% c("Land","Land-Ocean")))+
  geom_sf(aes(fill=Continent))+
  scale_fill_manual(values = colorvec)+
 geom_sf_text(aes(label = Name), size = 3, position = position_jitter(height = 3, width = 3, seed = 1))+
  #geom_sf_text(aes(label = Acronym), size = 2)+
  theme_void()+
  theme(legend.position = "bottom")

#ggsave("ipcc_regions.png", width = 12, height = 8, units = "cm", dpi = "retina")

3 Aridity Index

comp <- merge(select(era5vA, c("lon","lat","source","Continent","z","Rn","ET0","AI","cat.AI")), 
                     select(wdcvA, c("lon","lat","source","Continent","z","Rn","ET0","AI","cat.AI")), by = c("lon","lat")) %>%
        merge(select(cmip6vA, c("lon","lat","source","Continent","z","ET0","AI","cat.AI")), by = c("lon","lat"))

3.1 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000


calib <- rbind(cmip6vA, era5vA, wdcvA)

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"

ggarrange(
  
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "Wordclim, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "ERA5, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "CMIP6, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    scale_fill_manual(values = col.cat)+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "left")

)

3.2 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000 with country borders


calib <- rbind(cmip6vA, era5vA, wdcvA)

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"

ggarrange(
  
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey20")+
    labs(title = "Wordclim, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey20")+
    labs(title = "ERA5, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey20")+
    labs(title = "CMIP6, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    scale_fill_manual(values = col.cat)+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top"),

nrows = 4, ncol = 1

)


bpr <- c(-400, 0, 400, 800, 1200, 1600, 5000)
colvecrdbu <-  brewer.pal(n = 11, name ="RdYlBu")

col.pr <- c("(-400,0]" = colvecrdbu[5], "(0,400]"= colvecrdbu[6],"(400,800]"= colvecrdbu[7], "(800,1.2e+03]" = colvecrdbu[8], "(1.2e+03,1.6e+03]"= colvecrdbu[9], "(1.6e+03,5e+03]" = colvecrdbu[10])

calib$pr.cat <- calib$pr %>% cut(breaks = bpr) 

g.pr.wdc <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = pr.cat))+
    scale_fill_manual(values= col.pr, na.translate = F)+
    borders(colour= "grey20")+
    labs(title = "Annual precipitation in Wordclim", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.pr.era5 <- ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = pr.cat))+
      scale_fill_manual(values= col.pr, na.translate = F)+
    borders(colour= "grey20")+
    labs(title = "Annual precipitation in ERA5", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.pr.cmip6 <- ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = pr.cat))+
    scale_fill_manual(values= col.pr, na.translate = F)+
    borders(colour= "grey20")+
    labs(title = "Annual precipitation CMIP6", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

l.pr <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = pr.cat))+
    scale_fill_manual(values= col.pr, na.translate = F)+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top", legend.key.width = unit(2, "cm"))

#)
bpr <- c(-400, 0, 400, 800, 1200, 1600, 5000)
colvecrdbu <-  brewer.pal(n = 11, name ="RdYlBu")

col.et0 <- c("(-400,0]" = colvecrdbu[8], "(0,400]"= colvecrdbu[7],"(400,800]"= colvecrdbu[6], "(800,1.2e+03]" = colvecrdbu[5], "(1.2e+03,1.6e+03]"= colvecrdbu[4], "(1.6e+03,5e+03]" = colvecrdbu[3])

calib$et0.cat <- calib$ET0 %>% cut(breaks = bpr) 

g.et0.wdc <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    borders(colour= "grey20")+
    labs(title = "Annual evapotranspiration Wordclim", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.et0.era5 <- ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    borders(colour= "grey20")+
    labs(title = "Annual evapotranspiration ERA5", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.et0.cmip6 <- ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    borders(colour= "grey20")+
    labs(title = "Annual evapotranspiration CMIP6", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

l.et0 <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    geom_raster(aes(x = lon, y = lat), fill = "white")+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top", legend.key.width = unit(2, "cm"))
ggarrange(
  ggarrange(plotlist = list(g.pr.wdc,g.pr.era5, g.pr.cmip6, l.pr), nrow = 4, ncol = 1),
  ggarrange(plotlist = list(g.et0.wdc,g.et0.era5, g.et0.cmip6, l.et0), nrow = 4, ncol = 1),
  ncol = 2)

3.3 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000 with Antarctica


calib <- rbind(cmip6vA, era5vA, wdcvA)

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"

ggarrange(
  
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "Wordclim, 1970-2000", fill = "")+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "ERA5, 1970-2000", fill = "")+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "CMIP6, 1970-2000", fill = "")+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    scale_fill_manual(values = col.cat)+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "left")

)

calib$same.cat <- NA
calib$cat.ref <- NA

for(i in 1:nrow(calib)){
  loni = calib$lon[i]
  lati = calib$lat[i]
  catrefi <- subset(calib, lon == loni & lat == lati & source == "CMIP6")$cat.AI %>% as.character()
  cati <- calib$cat.AI[i] %>% as.character()
  calib$same.cat[i] <- ifelse(cati == catrefi, "y",
                              ifelse(cati!= catrefi, "n", 
                                     NA))
  calib$cat.ref[i] <- catrefi
}

write.table(calib, "calib.txt")
# table_calib_c <- calib %>% select(c("Continent","source","cat.AI")) %>%
#   reshape2::melt(id.vars = c("Continent","source")) %>%
#   group_by(Continent,source,value) %>%
#   summarise(count = n()) %>% mutate(perc = round(count/sum(count)*100,1)) %>%
#   ungroup()
# 
# table_calib_c$value[is.na(table_calib_c$value) == T] <- "NA"
# 
# ggplot(table_calib_c)+
#   geom_col(aes(x = Continent, y = perc, fill = value, position = "fill"))+
#   scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","NA"))+
#   theme_minimal()+
#   theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
#   facet_grid(rows = vars(source))+
#   labs(x = "Continent", y = "%", fill = "Aridity\ncategory", y = "Proportion of aridity category")
# 


table_calib <- calib %>% select(c("Continent","Name","Acronym","source","cat.AI")) %>%
  reshape2::melt(id.vars = c("Continent","Name","Acronym","source")) %>%
  filter(!is.na(value)) %>%
  group_by(Continent,Name,Acronym,source,value) %>%
  summarise(count = n()) %>% mutate(perc = round(count/sum(count)*100,1)) %>%
  ungroup() %>% filter(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC"))

table_calib$value[is.na(table_calib$value) == T] <- "NA"
table_calib$Continent[which(table_calib$Continent == "EUROPE-AFRICA")] <- "EUROPE"
table_calib$Continent[which(table_calib$Continent == "CENTRAL-AMERICA")] <- "CENTRAL\nAMERICA"
table_calib$Continent[which(table_calib$Continent == "NORTH-AMERICA")] <- "NORTH\nAMERICA"
table_calib$Continent[which(table_calib$Continent == "SOUTH-AMERICA")] <- "SOUTH\nAMERICA"

table_calib$Continent <- factor(table_calib$Continent, levels = c("AFRICA","EUROPE","ASIA","OCEANIA","NORTH\nAMERICA","CENTRAL\nAMERICA","SOUTH\nAMERICA"))

ggplot(table_calib)+
  geom_col(aes(x = Acronym, y = perc, fill = value, position = "fill"))+
  scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","NA"))+
  theme_minimal(base_size = 12)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), 
        legend.position = "bottom", 
        strip.text.x = element_text(size = 6))+
  facet_grid(source ~ Continent, switch = "y", scales = "free_x", space = "free")+
  labs(x = "", y = "%", fill = "Aridity\ncategory")


ggplot(table_calib)+
  geom_col(aes(x = Acronym, y = count, fill = value, position = "fill"))+
  scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","NA"))+
  scale_y_continuous(position = "right")+
  theme_minimal(base_size = 12)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), 
        legend.position = "bottom", 
        strip.text.x = element_text(size = 6))+
  facet_grid(source ~ Continent, switch = "y", scales = "free_x", space = "free")+
  labs(x = "", y = "Count", fill = "Aridity\ncategory")

library(ggrepel)

pie_calib <- table_calib %>% group_by(source, value) %>% summarise(gridcells = sum(count)) %>% mutate(percent = round(gridcells/sum(gridcells)*100, 1), lab.y = (rev(cumsum(rev(gridcells)))) - gridcells*0.5)

ggplot(pie_calib)+geom_bar(aes(x = "", y = gridcells, fill = value), width = 1, stat = "identity")+
  coord_polar("y", start = 0)+
  scale_fill_manual(values = col.cat, na.translate = T)+
  facet_grid(cols = vars(source))+
  geom_label_repel(aes(x = "", y = lab.y, label = percent), nudge_x = 0.5)+
  labs(fill = "")+
  theme_void()+theme(legend.position = "bottom")

4 Number of inhabitants by region

density <- raster("Population_density/gpw_v4_population_density_rev11_2020_1_deg.tif")
plot(density)

popcount <- raster("Population_density/gpw_v4_population_count_rev11_2020_1_deg.tif")
plot(popcount)


land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land

popcount.df <- popcount %>%
  projectRaster(to = land_mask) %>% 
  as.data.frame(xy = T) %>% setNames(c("lon","lat","pop")) %>% 
  merge(land_mask.df, by = c("lon","lat")) %>% filter(lm == 1)

era5vApop <- merge(era5vA, popcount.df, by = c("lon","lat")) 

tab.pop.by.cat <- era5vApop %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(cat.AI, Continent) %>%
  summarise(pop = sum(pop)) 

write.csv(tab.pop.by.cat, "population.count.era5.csv")
knitr::kable(tab.pop.by.cat)
cat.AI Continent pop
Arid AFRICA 9.407395e+07
Arid ASIA 2.608739e+08
Arid CENTRAL-AMERICA 1.102838e+07
Arid EUROPE 3.523270e+05
Arid EUROPE-AFRICA 3.354203e+07
Arid NORTH-AMERICA 8.275538e+05
Arid OCEANIA 2.483547e+05
Arid SOUTH-AMERICA 2.714742e+06
Cold ASIA 9.005723e+06
Cold EUROPE 6.822776e+06
Cold NORTH-AMERICA 8.555420e+05
Cold POLAR NA
Cold SOUTH-AMERICA 6.222267e+04
Dry subhumid AFRICA 1.131375e+08
Dry subhumid ASIA 3.724373e+08
Dry subhumid CENTRAL-AMERICA 2.847546e+07
Dry subhumid EUROPE 5.195230e+06
Dry subhumid EUROPE-AFRICA 5.933504e+07
Dry subhumid NORTH-AMERICA 9.392606e+06
Dry subhumid OCEANIA 7.027896e+05
Dry subhumid SOUTH-AMERICA 1.297947e+07
Humid AFRICA 6.854010e+08
Humid ASIA 3.055659e+09
Humid CENTRAL-AMERICA 1.462475e+08
Humid EUROPE 4.999673e+08
Humid EUROPE-AFRICA 1.452014e+08
Humid NORTH-AMERICA 2.915502e+08
Humid OCEANIA 1.529713e+07
Humid POLAR 4.329376e+02
Humid SOUTH-AMERICA 3.724919e+08
Hyperarid AFRICA 5.328568e+07
Hyperarid ASIA 1.970440e+07
Hyperarid EUROPE-AFRICA 4.146658e+07
Hyperarid SOUTH-AMERICA 5.088920e+05
Semi-arid AFRICA 2.154724e+08
Semi-arid ASIA 5.713649e+08
Semi-arid CENTRAL-AMERICA 3.973955e+07
Semi-arid EUROPE 1.841667e+06
Semi-arid EUROPE-AFRICA 8.619193e+07
Semi-arid NORTH-AMERICA 2.806716e+07
Semi-arid OCEANIA 3.403064e+06
Semi-arid SOUTH-AMERICA 2.872740e+07
NA AFRICA 3.321493e+04
NA ASIA 1.762974e+07
NA CENTRAL-AMERICA 7.132541e+05
NA EUROPE 1.481967e+05
NA EUROPE-AFRICA 1.876606e+06
NA NORTH-AMERICA 1.403549e+03
NA OCEANIA 3.896168e+04
NA POLAR NA
NA SOUTH-AMERICA 2.214226e+05

5 Comparison of datasets

vars <- c("lon","lat","Continent","Name","Acronym","tas","pr","sfcWind","Rn","ea","ET0","AI","cat.AI")
merge.vars <- c("lon","lat","Continent","Name","Acronym")
calib.wide <- select(era5vA, vars) %>%
               merge(select(wdcvA, vars), 
                     by = merge.vars) %>%
          merge(select(cmip6vA, vars), 
                by = merge.vars) %>%
          setNames(c(merge.vars,
                     "tas.era5","pr.era5", "sfcWind.era5","Rn.era5","ea.era5","ET0.era5","AI.era5","cat.AI.era5",
                     "tas.wdc","pr.wdc", "sfcWind.wdc","Rn.wdc","ea.wdc","ET0.wdc","AI.wdc","cat.AI.wdc",
                     "tas.cmip6","pr.cmip6", "sfcWind.cmip6","Rn.cmip6","ea.cmip6","ET0.cmip6","AI.cmip6","cat.AI.cmip6"))

5.1 Scatterplots with Antarctica

5.1.1 ERA5 compared to WDC

Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.

5.1.2 ERA5 and WDC

ggarrange(

ggplot(calib.wide,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

knitr::kable(data.frame(precipitation = cor(calib.wide$pr.wdc, calib.wide$pr.era5, use = "pairwise.complete.obs"),
                        temperature = cor(calib.wide$tas.wdc, calib.wide$tas.era5, use = "pairwise.complete.obs"),
                        wind = cor(calib.wide$sfcWind.wdc, calib.wide$sfcWind.era5, use = "pairwise.complete.obs"),
                        Rn = cor(calib.wide$Rn.wdc, calib.wide$Rn.era5, use = "pairwise.complete.obs"),
                        ea = cor(calib.wide$ea.wdc, calib.wide$ea.era5, use = "pairwise.complete.obs"),
                        ET0 = cor(calib.wide$ET0.wdc, calib.wide$ET0.era5, use = "pairwise.complete.obs")), 
             digits = 2, caption = "Correlation coefficients between ERA5 and Worldclim") %>%
  kableExtra::kable_styling(bootstrap_options = "bordered")
Correlation coefficients between ERA5 and Worldclim
precipitation temperature wind Rn ea ET0
0.89 0.99 0.68 0.67 1 0.98
ggplot(calib.wide) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none")


ggarrange(

ggplot(calib.wide) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.1.3 CMIP6 and WDC

ggarrange(

ggplot(calib.wide,aes(x = tas.wdc))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2
)

5.1.4 CMIP6 and ERA5


ggarrange(
  
ggplot(calib.wide,aes(x = tas.era5))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6", col = "Dataset")+
  scale_color_manual(values = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.era5))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6", col = "Dataset")+
  scale_color_manual(values = c(colorvec[8]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.era5))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6", col = "Dataset")+
  scale_color_manual(values = c(colorvec[8]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.era5))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
   stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, W/MJ/m2/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.era5))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5", y = "ea, kPa, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.era5))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2)

5.2 Scatterplots without Antarctica

calib.noant <- calib.wide %>% filter(lat > -60)

5.2.1 ERA5 compared to WDC

Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.

5.2.2 ERA5 and WDC

ggarrange(

ggplot(calib.noant,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

ggplot(calib.noant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none")


ggarrange(

ggplot(calib.noant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.2.3 CMIP6 and WDC

ggarrange(

ggplot(calib.noant,aes(x = tas.wdc))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.wdc))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.wdc))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.wdc))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.wdc))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.wdc))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2
)

5.2.4 CMIP6 and ERA5


ggarrange(
  
ggplot(calib.noant,aes(x = tas.era5))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6", col = "Dataset")+
  scale_color_manual(values = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.era5))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6", col = "Dataset")+
  scale_color_manual(values = c(colorvec[8]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.era5))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6", col = "Dataset")+
  scale_color_manual(values = c(colorvec[8]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.era5))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
   stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, W/MJ/m2/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.era5))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5", y = "ea, kPa, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.era5))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2)

5.3 Scatterplots only Antarctica

calib.ant <- calib.wide %>% filter(lat < -60)

5.3.1 ERA5 compared to WDC

Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.

5.3.2 ERA5 and WDC

ggarrange(

ggplot(calib.ant,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

ggplot(calib.ant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none")


ggarrange(

ggplot(calib.ant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.3.3 CMIP6 and WDC

ggarrange(

ggplot(calib.ant,aes(x = tas.wdc))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.wdc))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.wdc))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.wdc))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.wdc))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.wdc))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2
)

5.3.4 CMIP6 and ERA5


ggarrange(
  
ggplot(calib.ant,aes(x = tas.era5))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6", col = "Dataset")+
  scale_color_manual(values = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.era5))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6", col = "Dataset")+
  scale_color_manual(values = c(colorvec[8]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.era5))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6", col = "Dataset")+
  scale_color_manual(values = c(colorvec[8]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.era5))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
   stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, W/MJ/m2/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.era5))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5", y = "ea, kPa, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.era5))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2)

5.4 Aridity category difference

5.4.1 ERA5 and WDC

calib.wide$same.cat.wdc <- with(calib.wide, ifelse(cat.AI.era5 == cat.AI.wdc, "y","n"))
calib.wide$change.cat.wdc <- with(calib.wide, paste(cat.AI.era5, cat.AI.wdc, sep = " to "))
table(calib.wide$change.cat.wdc)
## 
##                 Arid to Arid                Arid to Humid 
##                         1359                            1 
##            Arid to Hyperarid            Arid to Semi-arid 
##                           18                          155 
##                 Cold to Arid                 Cold to Cold 
##                           19                        10101 
##         Cold to Dry subhumid                Cold to Humid 
##                            9                          290 
##            Cold to Semi-arid         Dry subhumid to Arid 
##                           20                           23 
##         Dry subhumid to Cold Dry subhumid to Dry subhumid 
##                            5                          205 
##        Dry subhumid to Humid    Dry subhumid to Semi-arid 
##                          209                          271 
##                Humid to Arid                Humid to Cold 
##                           30                          520 
##        Humid to Dry subhumid               Humid to Humid 
##                          347                         5436 
##           Humid to Semi-arid            Hyperarid to Arid 
##                          210                          230 
##            Hyperarid to Cold       Hyperarid to Hyperarid 
##                            2                          726 
##                   NA to Cold            Semi-arid to Arid 
##                          570                          181 
##            Semi-arid to Cold    Semi-arid to Dry subhumid 
##                            3                          173 
##           Semi-arid to Humid       Semi-arid to Semi-arid 
##                           71                         1077
calib.wide$dryer.wdc <- with(calib.wide, ifelse(change.cat.wdc == "Arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Arid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == " Arid to Dry subhumid", "dryer", 
                                ifelse(change.cat.wdc == "Arid to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Arid to Cold", "dryer",       
                                ifelse(change.cat.wdc == "Cold to Humid", "dryer",
                                ifelse(change.cat.wdc == "Cold to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Cold","dryer",
                                ifelse(change.cat.wdc == "Humid to Dry subhumid","wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Humid", "dryer",
                                NA))))))))))))))))))))))))

ggarrange(
ggplot(filter(calib.wide, same.cat.wdc == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.wdc))+
  scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
  borders()+
  labs(fill = "ERA5 compared to WDC", title = "Change of category")+
  ylim(-55,90)+
  theme_void(base_size = 15)+
  theme(legend.position = "bottom"),

ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat,  fill = AI.era5-AI.wdc))+
  binned_scale(aesthetics = "fill", breaks = c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), 
               palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")),
               guide = guide_legend(label.theme = element_text(angle = 0))
               )+
  labs(title = "Changes in AI", fill = "ERA5 compared to WDC")+
  theme_void(base_size = 15)+ylim(-55,90)+
  theme(legend.position = "bottom"),
ncol = 2)

ggplot() + 
  geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc < 277 & pr.era5 - pr.wdc > - 98), aes(x=lon, y = lat,  fill = pr.era5-pr.wdc))+
  geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc > 277), aes(x=lon, y = lat), fill = "#44FF44")+
   geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc < - 98), aes(x=lon, y = lat), fill = "#FF4444")+
  labs(title = "Most extreme differences in precipitation", subtitle = "In red, the 10% most extreme gricells in which Worldclim is wetter than ERA5; \nin green, the 10% most extreme gridcells in which ERA5 is wetter than Worldclim", fill = "precipitation in ERA5 - precipitations in WDC,\nmm/year")+
  theme_void(base_size = 15)+
  theme(legend.position = "bottom")

5.4.2 CMIP6 and WDC

calib.wide$same.cat.wdc <- with(calib.wide, ifelse(cat.AI.cmip6 == cat.AI.wdc, "y","n"))
calib.wide$change.cat.wdc <- with(calib.wide, paste(cat.AI.cmip6, cat.AI.wdc, sep = " to "))
table(calib.wide$change.cat.wdc)
## 
##                 Arid to Arid         Arid to Dry subhumid 
##                         1099                            9 
##            Arid to Hyperarid            Arid to Semi-arid 
##                           83                          150 
##                 Cold to Arid                 Cold to Cold 
##                           25                        10491 
##         Cold to Dry subhumid                Cold to Humid 
##                            5                          588 
##            Cold to Semi-arid         Dry subhumid to Arid 
##                           23                           36 
## Dry subhumid to Dry subhumid        Dry subhumid to Humid 
##                          147                          196 
##    Dry subhumid to Hyperarid    Dry subhumid to Semi-arid 
##                            1                          350 
##                Humid to Arid                Humid to Cold 
##                           77                          138 
##        Humid to Dry subhumid               Humid to Humid 
##                          427                         5112 
##           Humid to Hyperarid           Humid to Semi-arid 
##                            1                          435 
##            Hyperarid to Arid            Hyperarid to Cold 
##                           85                            2 
##       Hyperarid to Hyperarid                   NA to Cold 
##                          652                          570 
##            Semi-arid to Arid    Semi-arid to Dry subhumid 
##                          520                          146 
##           Semi-arid to Humid       Semi-arid to Hyperarid 
##                          111                            7 
##       Semi-arid to Semi-arid 
##                          775
calib.wide$dryer.wdc <- with(calib.wide, ifelse(change.cat.wdc == "Arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Arid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == " Arid to Dry subhumid", "dryer", 
                                ifelse(change.cat.wdc == "Arid to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Arid to Cold", "dryer",       
                                ifelse(change.cat.wdc == "Cold to Humid", "dryer",
                                ifelse(change.cat.wdc == "Cold to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Cold","dryer",
                                ifelse(change.cat.wdc == "Humid to Dry subhumid","wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Humid", "dryer",
                                NA))))))))))))))))))))))))

ggarrange(
ggplot(filter(calib.wide, same.cat.wdc == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.wdc))+
  scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
  borders()+
  labs(fill = "CMIP6 compared to WDC", "Change of category")+
  ylim(-55,90)+
  theme_void(base_size = 15)+
  theme(legend.position = "bottom"),

ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat,  fill = AI.cmip6-AI.wdc))+
   borders()+
  binned_scale(aesthetics = "fill", breaks = c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")), 
                                                                                                              guide = guide_legend(label.theme = element_text(angle = 0)))+
  labs(title = "Change in AI", fill = "CMIP6 compared to WDC")+
  theme_void(base_size = 15)+ylim(-55,90)+
  theme(legend.position = "bottom"),

ncol = 2)

5.5 CMIP6 and ERA5

calib.wide$same.cat.era5 <- with(calib.wide, ifelse(cat.AI.cmip6 == cat.AI.era5, "y","n"))
calib.wide$change.cat.era5 <- with(calib.wide, paste(cat.AI.cmip6, cat.AI.era5, sep = " to "))
table(calib.wide$change.cat.era5)
## 
##                 Arid to Arid         Arid to Dry subhumid 
##                          989                            6 
##                Arid to Humid            Arid to Hyperarid 
##                            2                          230 
##            Arid to Semi-arid                 Cold to Cold 
##                          114                        10377 
##         Cold to Dry subhumid                Cold to Humid 
##                           10                          745 
##         Dry subhumid to Arid Dry subhumid to Dry subhumid 
##                           14                          169 
##        Dry subhumid to Humid    Dry subhumid to Semi-arid 
##                          219                          328 
##                Humid to Arid                Humid to Cold 
##                            5                           62 
##        Humid to Dry subhumid               Humid to Humid 
##                          382                         5456 
##           Humid to Semi-arid            Hyperarid to Arid 
##                          285                           14 
##       Hyperarid to Hyperarid                     NA to NA 
##                          725                          570 
##            Semi-arid to Arid    Semi-arid to Dry subhumid 
##                          511                          146 
##           Semi-arid to Humid       Semi-arid to Hyperarid 
##                          121                            3 
##       Semi-arid to Semi-arid 
##                          778
calib.wide$dryer.era5 <- with(calib.wide, ifelse(change.cat.era5 == "Arid to Humid", "dryer",
                                ifelse(change.cat.era5 == "Arid to Semi-arid", "dryer",
                                ifelse(change.cat.era5 == " Arid to Dry subhumid", "dryer", 
                                ifelse(change.cat.era5 == "Arid to Hyperarid", "wetter",
                                ifelse(change.cat.era5 == "Cold to Humid", "dryer",
                                ifelse(change.cat.era5 == "Dry subhumid to Cold", "dryer",
                                ifelse(change.cat.era5 == "Dry subhumid to Humid", "dryer",
                                ifelse(change.cat.era5 == "Dry subhumid to Arid", "wetter",
                                ifelse(change.cat.era5 == "Dry subhumid to Semi-arid", "wetter",
                                ifelse(change.cat.era5 == "Humid to arid", "wetter",
                                ifelse(change.cat.era5 == "Humid to Semi-arid", "wetter",
                                ifelse(change.cat.era5 == "Humid to Cold","dryer",
                                ifelse(change.cat.era5 == "Humid to Dry subhumid","wetter",
                                ifelse(change.cat.era5 == "Semi-arid to Dry subhumid", "dryer",
                                ifelse(change.cat.era5 == "Semi-arid to Arid", "wetter",
                                ifelse(change.cat.era5 == "Semi-arid to Humid", "dryer",
                                ifelse(change.cat.era5 == "Hyperarid to Arid", "dryer",
                                ifelse(change.cat.era5 == "Hyperarid to Semi-arid", "dryer",
                                NA)))))))))))))))))))

ggarrange(
ggplot(filter(calib.wide, same.cat.era5 == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.era5))+
  scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
  borders()+
  labs(fill = "CMIP6 compared to ERA5", "Change of aridity category")+
  ylim(-55,90)+
  theme_void()+
  theme(legend.position = "bottom"),


ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat,  fill = AI.cmip6-AI.era5))+
   borders()+
  binned_scale(aesthetics = "fill", breaks =c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")))+
  labs(title = "Change of AI", fill = "CMIP6 compared to ERA5")+
  theme_void()+ylim(-55,90)+
  theme(legend.position = "bottom"),
ncol = 2)

6 NDVI

# Download Gimms NDVI
library(gimms)
library(gganimate)

land_mask <- raster("Masks/land_sea_mask_1degree.nc4")
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land

#ndvi.ref <- downloadGimms(x= 1981,y= 2000,dsn = "NDVI")

list.ndvi <- list.files(path = "NDVI", pattern = "ndvi.*\\.nc4$", full.names = T)[1:2]

ndvi.rast <- rasterizeGimms(list.ndvi, keep = 0) %>% mean(na.rm = T) %>% projectRaster(land_mask) # rasterize Gimms does the quality control and removes NAs
ndvi.df <- ndvi.rast %>% as.data.frame(xy = T) %>% setNames(c("lon","lat","ndvi")) %>% merge(land_mask.df, by = c("lon","lat")) %>% filter(lm != 0)

write.table(ndvi.df, "NDVI/ndvi.df.txt")
ndvi.df <- read.table("NDVI/ndvi.df.txt") 

ggplot(subset(ndvi.df, lm == 1))+geom_raster(aes(x=lon, y = lat, fill = ndvi))+
  scale_fill_continuous_sequential(palette = "Greens 3", rev = T)+
  ylim(-55,90)+labs(fill = "NDVI")+
  theme_void()+theme(legend.position = "bottom")


calibn <- merge(calib, ndvi.df, by = c("lon","lat"))

ggplot(subset(calibn, !is.na(cat.AI)))+geom_boxplot(aes(x=cat.AI, y = ndvi, col = cat.AI))+facet_grid(rows = vars(source))+
  scale_color_manual(values = col.cat)+
  labs(x = "", y = "NDVI")+
  theme_minimal()+theme(legend.position = "none")


aov(ndvi~cat.AI, data = group_by(calibn, source)) %>% summary()
##                Df Sum Sq Mean Sq F value Pr(>F)    
## cat.AI          5   1026  205.17    8708 <2e-16 ***
## Residuals   43200   1018    0.02                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 23577 observations effacées parce que manquantes

range.ndvi <- calibn %>% filter(!is.na(cat.AI)) %>% group_by(source, cat.AI) %>% 
  summarise(range.ndvi = paste(round(quantile(ndvi, na.rm = T)[c(2,4)],2), collapse = " to ")) %>%
  cast(source~cat.AI)

cor.ndvi <- calibn %>% filter(!is.na(cat.AI)) %>% group_by(source, cat.AI) %>% 
  summarise(pearson.r = round(cor(ndvi, AI, use = "pairwise.complete"),2)) %>%
  cast(source~cat.AI)

knitr::kable(range.ndvi, caption = "NDVI range by source and aridity category")
NDVI range by source and aridity category
source Cold Hyperarid Arid Semi-arid Dry subhumid Humid
CMIP6 0.3 to 0.47 0.09 to 0.11 0.1 to 0.19 0.17 to 0.36 0.24 to 0.49 0.4 to 0.69
ERA5 0.27 to 0.47 0.09 to 0.11 0.11 to 0.21 0.19 to 0.37 0.24 to 0.48 0.42 to 0.68
Worldclim 0.29 to 0.47 0.09 to 0.11 0.1 to 0.19 0.18 to 0.34 0.28 to 0.44 0.45 to 0.7
knitr::kable(cor.ndvi, caption = "Correlation between NDVI and AI by source and aridity category")
Correlation between NDVI and AI by source and aridity category
source Cold Hyperarid Arid Semi-arid Dry subhumid Humid
CMIP6 -0.08 0.15 0.36 0.38 -0.01 0.31
ERA5 -0.16 0.08 0.47 0.30 0.08 0.38
Worldclim 0.05 0.05 0.49 0.39 0.14 0.52
knitr::knit_exit()

6.0.1 Maps CMIP6 ERA5 // WDC

In red: when CMIP6 or ERA5 are warmer than Wordlclim

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.era5-tas.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: ERA5 - WDC")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.cmip6-tas.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: CMIP6 - WDC")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when CMIP6 or ERA5 are wetter than Worldclim

cowplot::plot_grid(
  
  ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.era5-pr.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: ERA5 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.cmip6-pr.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: CMIP6 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when wind speed is higher for CMIP6 or ERA5 compared to Worldclim

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.era5-sfcWind.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: ERA5 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.cmip6-sfcWind.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: CMIP6 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

6.0.2 Maps CMIP6 WDC // ERA5

In red: when CMIP6 or WDC are warmer than ERA5

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.wdc-tas.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: WDC - ERA5")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.cmip6-tas.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: CMIP6 - ERA5")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when CMIP6 or WDC are wetter than ERA5

cowplot::plot_grid(
  
  ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.wdc-pr.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: WDC - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.cmip6-pr.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: CMIP6 - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when wind speed is higher for CMIP6 or WDC compared to ERA5

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.wdc-sfcWind.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: WDC - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.cmip6-sfcWind.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: CMIP6 - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)